# -*- coding: utf-8 -*-
from sklearn.preprocessing import StandardScaler
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from scipy.stats import pearsonr

def partial_correlation_matrix(df):
    '''
    计算偏相关系数
    '''
    feature_num = df.shape[1]
    feature_name = df.columns
    partial_corr_matrix = np.zeros((feature_num, feature_num))
    # 计算偏相关矩阵
    for i in range(feature_num):
        x1 = df.iloc[:, i]
        for j in range(feature_num):
            if i == j:
                partial_corr_matrix[i, j] = 1
            elif j < i:
                partial_corr_matrix[i, j] = partial_corr_matrix[j, i]
            else:
                x2 = df.iloc[:, j]
                df_control = df.drop(columns=[feature_name[i], feature_name[j]], axis=1)
                L = LinearRegression().fit(df_control, x1)
                Lx = L.predict(df_control)
                x1_prime = x1 - Lx
                
                L = LinearRegression().fit(df_control, x2)
                Lx = L.predict(df_control)
                x2_prime = x2 - Lx
                partial_corr_matrix[i, j] = pearsonr(x1_prime, x2_prime)[0]
    return partial_corr_matrix

def prepare_data(time_series, situation):
    '''
    数据准备，导出时序数据
    '''
    # 总周数
    week_num = abs(time_series['接警日期'].iloc[-1] - time_series['接警日期'].iloc[0]).days//7
    # 找出所有的事故类别
    problem_classes = list(np.unique(time_series['事件类别']))
    accident_sites = list(np.unique(time_series['事件所在的区域']))
    data = pd.DataFrame()
    # 
    for accident_site in accident_sites:
        # 某地点的情况
        situation_site = situation.loc[situation['编号']==accident_site]
        # 面积
        site_area = situation_site['面积（km2）'].values
        site_area = np.float(site_area)
        # 事件类别分类
        time_series_site = time_series.loc[time_series['事件所在的区域']==accident_site]
        time_series_class = time_series_site.groupby(time_series_site['事件类别'])\
                            .count().iloc[: ,1]
        data[f'地点{accident_site}事件密度s'] = time_series_class/(week_num*site_area)
        
    data = data.fillna(0)
    return data    
    



if __name__ == '__main__':
    path = r'../附件/附件2：某地消防救援出警数据.xlsx'
    time_series = pd.read_excel(path)
    path = r'../附件/附件1：各区域的人口、面积.xlsx'
    situation = pd.read_excel(path)
    data = prepare_data(time_series, situation)
    
    # 保存数据
    data.to_excel(r'../附件/事件密度-事件类型.xlsx')
    
    print('各类型事件在不同区域的方差： \n', data.std(axis=1))
 
    matrix = partial_correlation_matrix(data.T)
    print('各类型事件在不同区域的偏相关矩阵如下：\n', matrix)
